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Abstract 

This paper considers the multi-parametric hnear complementarity problem (pLCP) with suf- 
ficient matrices. The main result is an algorithm to find a polyhedral decomposition of the set 
►^ , of feasible parameters and to construct a piecewise affine function that maps each feasible pa- 

QQ ' rameter to a solution of the associated LCP in such a way that the function is afhne over each 

cell of the decomposition. The algorithm is output-sensive in the sense that its time complexity 
^\J ' is polynomial in the size of the input and linear in the size of the output, when the problem is 

non-degcncratc. We give a lexicographic perturbation technique to resolve degeneracy as well. 
C*"- ' Unlike for the non-parametric case, the resolution turns out to be nontrivial, and in particular, 

^— ^ , it involves linear programming (LP) duality and multi-objective LP. 
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1 Introduction 

Given a real square matrix M and a vector q, solving a linear complementarity problem (LCP) 
consists of finding two nonnegative vectors w and z that satisfy the conditions 

w-Mz = q, w>Q, z>0, uFz = Q . (1.1) 

This simply stated and well-studied problem has far-reaching applications that have been well- 
documented in the literature. Rather than give a survey here, the interested reader is referred to 
the books [22, 6]. 

Several authors have studied the properties of various parametric versions of this problem 
(e.g. [17, 2, 6, 19, 7, 26, 8]), but unless there are restrictions placed on the particular parametric 
LCP (pLCP) considered, it is in general unrealistic to expect an efficient computational algorithm. 
We here study the class of pLCPs where the matrix M is sufficient^ and the right hand side (the 
vector q in (1.1)) is allowed to vary within a given affine subspace S. The goal is then to compute 



^Sufficient matrices are defined in Section 3. 



functions z(-) and w{-) that map from the affine subspace S" to a solution for pLCP (1.1) whenever 
one exists. 

This class of pLCP includes the important cases of linear and convex quadratic programs, where 
parameters appear linearly in the cost and the right hand side of the constraints [22]. In recent years, 
there has been a great deal of interest in the control community in parametric programming due to 
the fact that an important class of control algorithms for constrained linear systems, called model 
predictive controllers (MPC), can be posed as parametric linear or quadratic programs. The offline 
solution of these parametric problems results in an explicit representation of the optimal control 
action, which in some cases allows the controller to be implemented on systems with sampling rates 
of milli- and micro-seconds instead of the traditional seconds and minutes [25, 14, 4]. A similar 
setup results when computing optimal policies in a dynamic programming framework for partially 
observable Markov decision processes [18]. 

While parametric programming is widely used for sensitivity analysis, it is also applied in several 
other applications. In [15] it was shown that polyhedral projection can be reduced to parametric 
linear programming in polynomial time and of course such projections have uses ranging from the 
computation of invariant sets [5] and force closures [23] to program analysis [24] and theorem prov- 
ing [13]. Polyhedral vertex and facet enumeration can also be posed as projection problems, and 
hence solved with the proposed pLCP approach [12], which as discussed below results in an output 
sensitive algorithm in the non-degenerate case (although not the most efficient one for this purpose). 

In [2] it was shown that if M is a sufficient matrix and S satisfies certain general position'^ 
assumptions, then z{-) and w{-) are unique piecewise affine functions that are defined over a polyhe- 
dral partition of a convex set. There is, however, no known efficient method of testing this general 
position assumption a priori and in fact, it is often not satisfied even in the simplest case when the 
pLCP models a parametric linear program [16]. 

This paper extends the result of [2] by removing the restrictive and untestable general position 
assumption, allowing the algorithm to operate on any pLCP which is defined by a sufficient matrix 
and an affine subspace. This is achieved through a lexicographic perturbation technique, which 
has the effect of symbolically shifting the affine subspace an infinitesimally small amount and into 
general position. We first demonstrate that this perturbation always results in a problem that 
is in fact in general position and hence has the favorable uniqueness and partitioning properties 
discussed above. The challenge then becomes one of doing calculations in this perturbed space. The 
main optimization problem that arises as a result of the perturbation is a linear program that is 
polynomially parameterized by a positive variable e. The decision problem to be tackled is then the 
determination of the behavior of this parametric problem as the parameter e tends to zero. Section 5 
discusses how this problem can be converted into a multi-objective linear program, which can then 
be solved efficiently. The proposed technique should be applicable to other algorithms that rely on 
lexicographic perturbation to handle degeneracy. 

The resulting algorithm has the strong property that its complexity is polynomial in the size of 
the input (the matrix M) and linear in the size of the output (the number of pieces in the piecewise- 
affine functions w and z). For this reason, we call the algorithm 'output sensitive', although it should 
be noted that the complexity of the functions w and z can be exponential in the worst case and 
that this complexity result is for the lexicographically shifted affine subspace, which may be more 
complex than the unshifted case. 

The reminder of the paper is organized as follows. Section 2 gives some basic notations and 
a formal definition of the parametric LCP. Section 3 provides some useful properties of pLCPs on 
sufficient matrices. Section 4 then presents the proposed method with a general position assumption. 



^General position is defined in Section 3.2. 



and then this is relaxed in Section 5 where the lexicographic perturbation is introduced. Finally, 
Section 6 analyzes the complexity of the algorithm. 

2 Parametric LCP, critical regions and their adjacency 

Let us first fix some useful notations for matrices. For a matrix A € M'"^" and a column index j € 
{1, . . . , n}, A.j G M"^ denotes the j-th column vector of A. Similarly, for a row index i G {1, . . . , m} 
Ai G M^^" denotes the i-ih. row vector of A. For a subset J C {1, . . . ,n}, A.j G M™^"^ denotes the 
matrix formed by the columns of A indexed by J, and for a vector v G M", vj denotes the vector 
formed by the components of v indexed by J. For / C {1, . . . ,m}, we denote by Ai G M the 
matrix formed by the rows of A indexed by /. 

Given a real square matrix M and a vector q of size n, the linear complementarity problem (LCP) 
is to find two nonnegative vectors w and z that satisfy 

w-Mz = q, w>Q, z > 0, uFz = . (2.2) 

In this paper, we consider the LCP (2.2) where the right-hand side is allowed to vary within some 
affine subspace. Specifically, the goal is to find two functions w{-) and z{-) that solve (2.2) over a 
given affine subspace S. 

Definition 2.1. Let Q G M"^"^ be a matrix of rank d, q ^ M" a vector and M G M"^" a matrix of 
order n. The functions w{-) and z{-) are a solution to the pLCP (2.3) if for every 9 G 0/, w^O) and 
z{9) satisfy the relations 

w{e) - Mz{e) =q + Q9 , (2.3a) 

w{9),z(9) >0 , (2.3b) 

w{9fz{9) = , (2.3c) 

where Qf C W^ is the set of feasible parameters 9, that is, those for which a solution to (2.3) exists. 

For the remainder of the paper we assume that the problem data M, q and Q are given and 
we define A G M"^^" to be the matrix [/ — M] . Consider the following system of linear equality 
constraints in non-negative variables 

Ax = q , X > . (2.4) 

A basis is a set B C {1, 2, . . . , 2n} such that \B\ = n and rank{A.B) = n; N := {1, . . . , 2n}\B is its 
complement and we call xb and xn the basic and non-basic variables respectively. Every basis B 
defines a basic solution to the linear system (2.4) 

XB = A~j^q , XN = . (2.5) 

A basis B is called complementary if |{z,z -|- n} n i?| = 1 for all i = 1, . . . ,n, and feasible if the 
associated basic solution satisfies the nonnegativity constraint in (2.4), i.e. A~^q > 0. Every 
complementary feasible basis defines a solution of the LCP (2.2), by setting {w , z ) = x . In the 
parametric case, each basis is feasible for a set of parameters, which leads to the notion of a critical 
region. 

Definition 2.2. The critical region TZb of a complementary basis B is defined as the set of all 
parameter values for which B is feasible, i.e., 

7^B := {9eM.'^\ Ariiq + Q9) > 0} . (2.6) 

A complementary basis B is called feasible for the pLCP (2.3) if TZb is nonempty. 



By definition critical regions are convex polyhedra contained in the set of feasible parameters 
Qf. Each feasible complementary basis B defines a solution of the pLCP for each 9 S TZb as 
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= 0, which is an affine function in IZb- As a result, if 0j can 

N 

be partitioned into a set of critical regions whose interiors are disjoint, then we have immediately a 
piecewise affine solution of pLCP (2.3) defined over these critical regions. 

In this paper we define a set of conditions under which such a partitioning can be achieved and 
introduce an efficient algorithm for this class of problems. The algorithm is based on the tracing 
of a graph whose nodes are the full-dimensional critical regions and whose edges are the pairs of 
adjacent regions (having a (d — 1) -dimensional intersection). 

Definition 2.3. Two critical regions TZi,TZ2 are called adjacent if their intersection T^i n 7^2 is of 
dimension d — 1. 

Definition 2.4. Let V be the set of complementary bases B of pLCP (2.3) such that TZb is full- 
dimensional and let E be the set of pairs of bases in V whose critical regions are adjacent. The 
graph G := (V, E) is called the critical region graph of the pLCP (2.3). 

The proposed algorithm enumerates all full-dimensional critical regions by tracing the above 
graph. This tracing requires that we are able to enumerate all neighbors of a given complementary 
basis. The following section discusses the properties of this graph and investigates restrictions on 
matrices M under which the neighbor search can be done efficiently. 

3 Well behaving matrix classes for parametric LCPs 

The goal of solving a parametric LCP is to compute functions w{-) and z{-) that satisfy (2.3) for all 
feasible values of the parameter 9 S 0/. As discussed in the introduction, linear complementarity 
problems include a very large set of difficult optimization problems and so we cannot hope for a 
solution in the general case. In this section, we identify classes of LCPs that are 'well-behaving', or 
that have properties which guarantee that the algorithm given in Section 4 will find a solution. 

The two key properties that will be needed are convexity of the feasible set 0/ and the existence 
of a "canonical" single- valued mapping from parameters to critical regions. The latter essentially 
means that the relative interiors of critical regions do not intersect. In this section we will formalize 
these notions and discuss a well-known matrix class that has the appropriate properties when the 
affine subspace S C M" of all possible right hand sides is the whole space M". In Section 3.2 we will 
then generalize this and give conditions such that these properties still hold when the right hand 
side is restricted to lie in some lower-dimensional affine subspace.^ 

3.1 Complementary cones 

We begin by describing the set of right hand sides q in (2.3) that are feasible for a given set of active 
constraints. 

For any index i G {!,..., 2n} we denote with i the complementary index of i, i.e. i = {i + n) 
mod 2n. For a set / C {!,..., 2n}, the set I is defined as the set of all complementary indices of 
elements in /. A set J C {!,..., 2n} is called complementary if i € J implies i ^ J . 



^Throughout the paper, we use the same notation regarding matrix classes as in [6] and we use the properties of 
each class proved there. At the end of the paper we append an auxiliary section, where the relevant definitions and 
theorems are mentioned. 



Definition 3.1. For any complementary set J, the cone C{J) := cone(A.j) is called a complementary 
cone (relative to M), where cone(r) denotes the cone of all nonnegative combinations of the columns 
of a matrix T. 

If i? is a complementary basis, then the complementary cone C{B) is full-dimensional, and 
conversely if the complementary cone C{J) is full-dimensional then the submatrix A.j has full rank, 
i.e. J is a complementary basis. For a complementary basis i?, we have 

C{B) = {y € M" I J^j^y > 0} . (3.7) 

In the remainder of the paper we will denote by /3 the matrix A"^ , where B is the considered basis. 
Therefore we wiU write C{B) = {y €W\f3y> 0}. 

One can see that for a given basis B, the cone C{B) is the set of all right hand sides that are 
feasible for LCP (2.2). We are interested in LCPs that have complementary cones with disjoint 
interiors and so we introduce the class of sufficient matrices, which has this property. 

Definition 3.2. A matrix M G M"^" is called column sufficient if it satisfies the implication 

[zi{Mz)i < for ah i] =^ [zi{Mz)i = for all i] . (3.8) 

The matrix M is called row sufficient if its transpose is column sufficient. If M is both column and 
row sufficient, then it is called sufficient. 

Remark 3.3. We note that both positive semidefinite (abbreviated by PSD) and P-matrices are 
sufficient. For a given matrix M it is possible to test in finite time whether it is sufficient, although 
no polynomial time test is currently known. 

The class of LCPs with sufficient matrices has been studied extensively, partly because this class 
appears to capture all critical structures for LCPs to behave nicely. In particular, this class admits 
many fruitful results ranging from combinatorial algorithms and duality [11, 10] to the efficient 
solvability by interior-point methods [20] . We will see that this class is ideal also for the investigation 
of parametric LCPs. We start with a key fact. 

Proposition 3.4 ([6, Theorem 6.6.6]). If M is a sufficient matrix, then the relative interiors of any 
two distinct complementary cones are disjoint. 

The union of all complementary cones forms a set known as the complementary range K[M). 
The complementary range is equal to the set of all right hand sides of the LCP for which a feasible 
solution exists [6] 

K{M) := {q \ the LCP (2.2) with matrix M and right hand side q is feasible} . (3.9) 

Proposition 3.5. If M is a sufficient matrix, then the complementary range K{M) is a convex 
polyhedral cone K{M) = cone([/ — M]). 

Proof. The statement follows from the fact that sufficient matrices are in Qo, see Theorem A. 10. D 

Remark 3.6. Throughout the paper we will draw upon the properties of two matrix classes exten- 
sively. The first class is the Qo-matrices, whose complementary range K[M) is a convex cone and 
the second is the fully semi-monotone matrices, denoted by Eq which have complementary cones 
that are all disjoint in their interiors. The class of sufficient matrices is contained in Qo n Eq and 
is perhaps the largest known subclass defined by a simple set of conditions, which is why sufficiency 
is assumed for the majority of the results in this paper. It should be noted, however, that many of 
the results hold under slightly relaxed assumptions. 



We will study now the adjacency relationship of complementary cones for the case of sufficient 
matrices. Specifically, since our goal is to compute the critical region graph Q, finding all neighbors 
of any given region is a crucial issue. We first look at the neighbors of a complementary cone that 
determine possible candidates for the neighbors for a critical region. 

Definition 3.7. Two complementary bases Bi and B2 are called adjacent if their cones C{Bi) and 
C{B2) are adjacent, that is, the dimension of C{Bi) n C{B2) is n — 1. 

The following lemma is important in narrowing down the candidates of the neighbor search. 

Lemma 3.8. //M € M"^" is a sufficient matrix and Bi and B2 are adjacent complementary bases, 
then \Bi nB2\>n-2. 

Proof. By the definition of adjacency the intersection C{Bi)r]C{B2) has dimension n—1 and therefore 
there exists a (7 G K{M) that lies in the relative interior of a facet of both complementary cones, which 
means that both basic solutions (^1,2:1), {w2,Z2) have exactly n — 1 strictly positive components. 
Recall that basic solutions can be stated as: 



Zl 



Bi 



W2 

Z2 



B'2 



Let J be a subset of Bi such that {wf,zf)j > and \J\ = n — 1. By Theorem A. 11, we have 
(^2 5-22 )j = 0. Since {w2 , Z2)b2 has exactly one zero component, at least n — 2 elements of J are 
not in B2 and therefore their complements are. This shows |i?i fl -B2J > n — 2. D 

Remark 3.9. It can be shown for P-matrices that two bases are adjacent if and only if they 
differ by exactly one element. This implies that the set of all complementary cones for a P-matrix 
LCP together with their faces forms a polyhedral complex. Unfortunately, this polyhedral complex 
property is not satisfied in general for sufficient matrices, nor in fact for the proper subclass of PSD 
matrices. More precisely, the intersection of two critical regions may not be a common face, see [17]. 

Lemma 3.10. Let M € R"^" be sufficient and B be a complementary basis. If B' = B\{i} U {i} 
is a basis then C{B) and C{B') intersect in their common facet C{B\{i}). Moreover no other full- 
dimensional complementary cones intersect the relative interior ofC{B\{i}), i.e. C{B') is the unique 
complementary cone adjacent to C{B) along this facet. 

Proof. Since B\{i} is a subset of B and B\ C{B\{i}) is a common facet of C{B) and C{B'). The 
second statement follows directly from the fact that the interior of any other complementary cone 
can intersect neither C{B) nor C{B'), since M is sufficient. D 

2 

Remark 3.11. Lemmas 3.8 and 3.10 imply that a complementary basis has at most n + " ^" 
adjacent complementary bases. 

Given a complementary basis B one can see that replacing any index i G B with its complement 
i preserves complementarity, i.e. B\{i} U {i} is still a complementary set. This operation is called 
a diagonal pivot. If we substitute two different indices i,j G B with their complements, then the 
operation is called an exchange pivot. Lemma 3.8 ensures that for a given basis B we can reach all 
adjacent bases by a single diagonal pivot or by a single exchange pivot operation. However, for some 
i £ B the set B U {i}\{i} may not be a basis, or for some pair {i,j) G B the basis B U {i,j}\{i,j} 
may not be adjacent to B. Therefore, in order to determine whether a set given by a diagonal or 
an exchange pivot is in fact an adjacent feasible basis we need a further condition. Such a condition 
can be easily derived from the dictionary of the basis B. 

6 



Definition 3.12. Given a complementary basis B and its complement A^ the matrix D := —A^ A.^ € 
M is called the dictionary of B. 

We begin by examining the diagonal pivot, for which a well-known adjacency condition can be 
derived. 

Fact 3.13. If B is a complementary basis, then for any i ^ B, the set B\{i'\ U {i} is a basis if and 
only if Dq / 0. 

We now consider the exchange pivot and derive necessary and sufficient conditions for adjacency, 
which are again based on examining elements of the dictionary. 

Proposition 3.14. Let M G M"^" be a sufficient matrix, B be a complementary basis and D be 
its dictionary. Consider the complementary basis B' = B\{i,j} U {i,j}, where i,j (z B are distinct. 
The following condition holds: 

dim{C{B\{i}) n C{B')) = n-\ -^^ D- = and Djj < 0. (3.10) 

Proof. Define a^ := —Dj^i, then the following holds: 

keB\{i} 

Let j G B\{i}, since aj / (we have assumed B' to be a basis) we can rewrite (3.11) as 



A.J = —iA.j- Yl '^kA.k ■ (3.12) 

"^' \ keB\{i,j} j 

Let us consider 

q{\) = Y ^kA.k , (3.13) 

keB\{i} 

which lies in the relative interior of C{B\{i}) if and only if A^ > for all k. 
We can express q{X) in following way by substituting (3.12) in (3.13): 

q{X)= Yj AfcAfc + Aj(l/ajAi- Y (^k/ajA.k) (3.14) 

fcGB\{i,i} fceB\{i,i} 

= Yj {>^k-ak/aj)A.k + Xj/ajA.-i (3.15) 

keB\{i,j} 

Sufficiency: if aj > then there exists a q{X) that lies in the relative interior of both facets C{B\{i}) 

andC(i?'\{i}). 

Necessity: since B' is a basis the unique way to express q{X) as a linear combination of the vector 

indexed by B' is (3.15). If aj < any q{X) E relint(C(-B\{i})) can not lie in C{B'). The case aj = 

is impossible since we have assumed B' to be a basis. D 

Corollary 3.15 follows directly from the proposition above and allows the detection of the bound- 
aries of the complementary range. 

Corollary 3.15. Let M G M"^" be sufficient, B be a complementary basis and denote A~^ as (5. 
Consider the facet C{B\{i]), for any i G B. The hyperplane aff(C(i?\{i})) = {y G W^\l3iy = 0} 
defines a facet of the complementary range K[M) if and only if D^j = and Di > 0. 

Remark 3.16. Lemma 3.8 and Proposition 3.14 are valid also for column sufficient matrices (i.e. 
it is not necessary that M is row sufficient). Lemma 3.10 holds for all fully semimonotone matrices. 
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3.2 Critical Domains 

We study now the parametric case where the right hand side of LCP (2.2) is restricted to he within 
some afhne subspace S := {q + QO \ 9 G W^}. We will see, under some assumptions on S, that the 
properties of the complementary cones discussed in the previous section still hold in this case. 

Definition 3.17. If i? is a complementary basis, then the critical domain Sb is the intersection of 
the affine subspace S with the complementary cone C{B) 

Sb := C{B) nS = {y\ AT^y > 0, y = Q^ + g, G M"^} . (3.16) 

Since we have assumed Q to be full column rank, the parametrisation q + Q9 is an invertible 
function and it is not hard to see that a critical domain is therefore the image of a critical region, i.e. 
Sb = QTIb + (1- Since the parametrisation is a bijection, Sb and TZb have the same combinatorial 
structure for any complementary basis B. In particular, we have: 

Remark 3.18. The inequality (3iy > is redundant in Sb if and only if f3i{Q0 + q) > is redundant 
in TZb, where /? = A~^ . 

We now define a key assumption, which will allow the extension of the properties of complemen- 
tary cones to critical domains. 

Definition 3.19. The afhne subspace S is said to lie in general position if for every complementary 
basis B the following condition holds 

S intersects C{B) =^ S intersects mt{C{B)) . (3-17) 

If a critical domain has dimension d = dim(S'), we simply say that it is full- dimensional. By the 
definition above, we have: 

Remark 3.20. If S lies in a general position, then every critical domain is either full-dimensional 
or empty. 

Proposition 3.21. // Af is sufficient and S lies in general position, then the relative interiors of 
critical domains Sbi and Sb2 ^-re disjoint for any two distinct complementary bases Bi and B2. 

Proof. The statement is a direct consequence of Proposition 3.4, i.e. int(C(i?i)) and int(C(i?2)) are 
disjoint, and from relint(S'Bj Q int(C(i?j)) for i = 1,2. D 

We denote the set of the feasible points of S by Sj. By (3.9) it follows that Sf = S (1 K{M). 

Corollary 3.22. If M G M"^" is sufficient then Sf is a convex polyhedron. 

Proof. The statement is a direct consequence of Proposition 3.5. D 

If M is a sufficient matrix and S is in general position, then Proposition 3.21 and Corollary 3.22 
ensure that the set K of nonempty critical domains defines a polyhedral decomposition of Sf in the 
sense that 

• each member P oi K is a convex polyhedron, 

• Up^kP = Sf, 

• dimP = d for all P ^ K, and 
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• dim(P n P') < d —1 for any two distinct members P and P' of K. 

It is important to note that the set K may not induce a polyhedral complex, i.e. the intersection of 
two critical domains may not be a common face. Nevertheless, because Sf is convex, we can define 
a graph structure of the decomposition which is connected. 

Definition 3.23. Let V be the set of complementary bases B of pLCP (2.3) such that Sb is full- 
dimensional and E consist of edges connecting each pair of bases in V whose critical domains are 
adjacent. The graph Q := (y,E) is called the critical domain graph of the pLCP (2.3). 

As stated above, each critical domain Sb is the image of a critical region TZb under the affine 
map 9 >-^ q + Q9, and a similar statement can be make for the feasible sets Sf = q + QQf. Since 
for each complementary basis B the critical domain Sb and the critical region TZb have the same 
combinatorial structure, the critical domain graph Q = {V, E) also defines the graph of critical regions 
and vice versa. In the discussion of the algorithm we will mostly consider only critical domains. 

Corollary 3.24. If M is a sufficient matrix and S lies in general position, then the graph of critical 
domains Q is connected. 

Proof. The statement follows directly from the convexity of Sf = K{M) f] S, which implies that 
between every pair of critical domains there exists a path in the graph of critical domains. D 

In the previous section we have seen that for the case of sufficient matrices we can reach all 
adjacent cones from any complementary cone with a single diagonal or a single exchange pivot 
operation. Assuming general position of S, this useful property also holds for critical domains. 

Proposition 3.25. If M is sufficient and S lies in general position, then for any two complementary 
bases Bi and B2 that have nonempty critical domains the following holds: If Sb^ ^.''^d Sb2 '^'^e adjacent 
in S then C{Bi) and C(i?2) o,re adjacent cones. 

Proof. If Sbj^ and Sb2 are adjacent critical domains, then their intersection is contained in C{Bi) n 
C{B2). If C(i?i) and C{B2) are not adjacent cones then there exists a complementary cone C adjacent 
to C{Bi) that contains C{Bi)i^C{B2). In this case Cr\S would be adjacent to Sbi and would overlap 
with the relative interior of Sb2^ which is a contradiction of Proposition 3.21. D 

Given a complementary feasible basis B, Proposition 3.25 ensures that all the critical domains 
adjacent to 5^ can be reached by exploring complementary bases adjacent to B. 

4 Description of the generic algorithm 

Now we are able to present an algorithm that enumerates all complementary bases whose critical 
domains define a polyhedral partition of Sf with the following two sets of assumptions: 

Assumption 4.1 (Regularity). The matrix M is sufficient and the matrix of the parametrisation 
Q G M"'^ has full column rank. 



Assumption 4.2 (General Position). The affine subspace S := {q + Q9\9 G M } lies in general 
position with respect to the complementary cones relative to M. 

Assumption 4.1 is essential for our algorithm to work, whereas Assumption 4.2 will be relaxed 
in the next section where an extension of the algorithm simulating general position for any given 
affine subspace via a symbolic perturbation is presented. 
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The proposed algorithm given in Algorithm 1 is based on a standard graph search procedure. It 
assumes a given function neighbors{B) , which returns all bases whose critical domains are adjacent 
to that of a given basis B. The validity follows immediately from the connectivity of the critical 
domain graph, Corollary 3.24. As input it takes a matrix M and an affine subspace S that satisfy 
the above assumptions, as well as an initial feasible complementary basis Bq such that Sbo is full- 
dimensional. The basis Bq is flagged as "unexplored" and added to the set of discovered bases B. 
In each iteration of the algorithm an unexplored basis B is selected from B, marked as "explored" 
and all bases that have adjacent critical domains are enumerated and added to B, marking the new 
bases as "unexplored". Once all bases in B have been explored, then we have found all bases with 
full-dimensional critical domains. 

Algorithm 1 Enumerate all critical domains by graph search 

Input: A feasible basis Bq with dim{SBo) = d, a sufficient matrix M E M"^'" and an affine 

subspace S that lies in general position. 
Output: The critical domain graph G = (B, E). 



7: 



Initialise the set of nodes B := {Bq} and edges E := 0. 
Flag Bq as "unexplored" . 

while there exists an unexplored basis B in B do 
Flag B as "explored" 

Bncw '■= neighbors{B) / /B^ew '■= neighbors'^ (B) if S does not lie in gen- 

eral position 
6: Flag each B' G Bncw\B as "unexplored" 
B:=BU i3ncw 

E :=EU {B, B') for each B' € ^ncw 
9: end w^hile 
10: return Q = {B, E) 

The remainder of this section describes how the results of the previous sections can be exploited to 
efficiently enumerate all adjacent critical domains of a given basis, i.e., how the function neighbors{-) 
can be properly implemented. The following section will then detail how the method can be extended 
so that the general position assumption can be relaxed. 

4.1 Neighborhood computation of a critical domain 

This section details a computational method that enumerates all bases that define adjacent crit- 
ical domains of a given basis, i.e. how the basis is "explored", under both Assumption 4.1 and 
Assumption 4.2. 

The function neighbors is given as Algorithm 2. Let S be a basis whose critical domain is full- 
dimensional. By Proposition 3.25, each adjacent critical domain must have a (d — 1) dimensional 
intersection with a facet of Sb- We begin therefore by first computing all facets of Sb and then by 
determining the critical domains that intersect each one. 

Given a complementary feasible basis B we determine which facets of C{B) define the facets of 
Sb by removing the redundant inequalities of Sb = {y ^ ^^ \ Py > 0, y = Qd + q}-, where /? := A~^ . 
The hyperplane hi = {y\f3iy = 0},i G B intersected with 5^ is a facet of Sb if there exists a 
y* G Sb such that f3jy* > for all j G B\{i} and (3iy* = 0. This fact relies on the general position 
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t* = max t 

s.t. -^jQe + t < pjq, Vj G B\{i} 
-PiQe = Piq 



assumption, Assumption 4.2. Therefore hi f] Sb is a facet of Sb if and only if the following LP: 

(4.18) 
PiQ9 = PiQ 

has an optimal value t* > strictly positive. 

Algorithm 2 Function neighbor s{B): returns all bases whose critical domains are adjacent to Sb- 
Input: A complementary basis i?, the matrix M and the affine subspace S. M is assumed to be 

sufficient and 5 to lie in general position. 
Output: The set B of complementary bases, whose critical domains are adjacent to Sb- 



1 

2 

3 

4 

5 

6 

7: 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 



B: 

/?: 
D 



A 



1 

■B 

A^ 



?xAr 



.B^N e.. 

for each i ^ B do 

if PiV > is non-redundant in Sb then 
if Dn > then 

add B\{i} U {i} to B 
else 

for each j £ B with Df^ < do 
B":=B\{i,j}Vj{iJ} 
if dim{SB n Sb") = d — \ then 

add B" to B 
end if 
end for 
end if 
end if 
end for 
return B 



//where iV = {l,...,2n}\S 
//by solving LP (4.18) 



//by solving LP (4.19) 



By solving LP (4.18) for each i € -B we can determine if C{B\{i}) defines a facet of Sb or not; 
see Line 5 of Algorithm 2. If it does, then the goal is to determine which bases, if any, have critical 
domains that intersect this facet. From the previous section, we saw that there are three possible 
cases: 

Diagonal pivot. If Df^ > then the cone C{B') defined by B' := B\{i\ U {%} is the unique 
complementary cone adjacent to C{B) along the facet C{B) n /ij, due to Lemma 3.10, and therefore 
Sb' is the unique critical domain adjacent to Sb along the facet Sb rihi. Since Sb' is nonempty (as 
Sb n hi is included in Sb') it is full-dimensional by Assumption 4.2. 

Boundary of K{M). C{B\{i}) is a facet of the complementary range and therefore no other 
complementary cones intersect it. From Corollary 3.15, this is the case when Dij = and Di > 0. 

Exchange pivot. By looking at the dictionary D oi B all complementary cones adjacent to C{B) 
that contain the index i can be determined (see Theorem 3.14). However, not all such cones intersect 
Sb with dimension d — 1 and for this reason we need to test for each cone C{B") adjacent to C{B) 
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(a) Two-dimensional slice of a three-dimensional ex- (b) Example of diagonal pivot: no adjacency check is 

ample. The cones C2 and C3 are both adjacent to needed. 

Ci along the same facet. However the afRne space S 
does not intersect C2 and therefore S fl C2 is not ad- 
jacent to 5 n Ci. This situation arises in the case of 
an exchange pivot and requires that adjacency much 
be checked by solving an LP. 

Figure 1: Adjacency of critical domains. See Example 1 for details. 

whether Sb" is adjacent to Sb, i-e. whether the condition dmi{SB'^SB") = d— 1 holds. Now assume 
that the basis B" := B\{i,j} U {i,j}, where j € B\{i}, defines such an adjacent complementary 
cone according to Theorem 3.14. In order to determine whether 5^" is adjacent to Sb we can solve 
following LP: 



(4.19) 



max 


t 




s.t. 


-PkQ0 + t < Pkq 


yk e B\{i} 




-l3'lQe + t < (3kq 


Vfc G B"\{j} 




-PiQO = p^q 






-P'JQ0 = P'Jq , 





where /3 






and (5" = A, 



As in the simpler case (4.18) above, the intersection Sb H Sb" 
has dimension d — 1 if and only if the optimal value of (4.19) is strictly positive. In this case 
Sb" = C{B") n 5 is nonempty and by the general position assumption, it is also full-dimensional. 



Remark 4.3. Since C{B) and C{B") are adjacent cones and the hyperplanes {y\f3iy = 0} and 
{ylP'j'y = 0} defines their shared facet, the two hyperplanes must be equivalent. Therefore in (4.19) 
one of the equality constraints —(3''Q6 = P'-q or —(3iQ9 = fiiq can be removed. 

Example 1. In Figure 1(a) a two-dimensional slice of three, three-dimensional cones is shown. The 
cones C2 and C3 are both adjacent to Ci along the same facet. However the affine subspace S does 
not intersect C2 and therefore 5 Pi C2 is not adjacent to S CiCi. 

In Figure 1(b) we consider the complementary basis Bi = {1, 2} and the cone C{Bi) = cone(/) = 
{y > 0}. The goal is to find the critical domain that is adjacent to Sbi = C{Bi) n S. The inequality 
yi > is not redundant in Sbi and therefore the hyperplane hi = {j/i = 0} defines a facet of Sbi • 
Since cone(— Mi, I2) n hi is equal to C{Bi) Dhi, we have that cone(— Mi, /2) CiS is adjacent to Sbi ■ 
The other inequality y2 >0 is redundant in Sb^ and therefore the critical domain cone(Ii, —M2) Pi S 
is not adjacent to Sbi ■ 
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5 Extension of the algorithm for S not in general position 

The previous section presented an algorithm that enumerates all feasible bases and returns the 
graph of critical domains. The algorithm works only under the assumption that the image of the 
parametrisation S lies in general position; Assumption 4.2. However, this assumption is not realistic 
and it is highly desirable to remove it. 

In the case of degeneracy (i.e. S is not in general position), Propositions 3.21 and 3.25 are no 
longer valid, as can be seen in Example 2. Therefore, during neighborhood computation it is not 
sufficient to explore only the adjacent complementary cones. In order to extend the algorithm to 
the degenerate case, we apply a symbolic perturbation technique (the lexicographic perturbation) 
which has the effect of shifting S into general position. 

The next subsection will demonstrate how to handle the perturbation for neighborhood compu- 
tation, in particular lines 5 and 11 of Algorithm 2. By using this technique we obtain a graph of 
critical domains G'' relative to the perturbed affine subspace 5"^, which can differ from the graph of 
critical domains Q relative to S. In particular, some full-dimensional critical domains in S^ may be 
non full-dimensional in S. We will see that there exists a subgraph Q of G^ that is a graph of critical 
domains relative to S and which can be obtained by postprocessing G^- 

Example 2. This example demonstrates the effect when a parametric LCP is not in general position. 
Consider the parametric LCP defined by the matrices 



M 



1 -1 
1 1 



, Q 



1 
-1 



and q 



A figure depicting the complementary cones relative to M and of the affine subspace S is shown in 
Figure 5. Let Bi = {1,2}, B2 = {2,3}, S3 = {3,4} and B4 = {1,4}. For notational simplicity, 
we denote by Ci = C{Bi) the complementary cones and by Si = S H Ci for i = 1, . . . ,4 the critical 
domains. Clearly S does not lie in general position because it intersects Ci,C3 and C4 on their 
boundary but not in their interiors. Proposition 3.21 is violated because Si is neither empty nor 
full- dimensional, and furthermore S3 and S4 are equal and hence relint(53) = relint(S'4). Theorem 
3.25 is violated because S2 and S4 are adjacent, but C2 and C4 are not. 



5.1 Lexicographic perturbation 

This section presents a well-known method that permits the perturbation of the image of the 
parametrisation S into general position and which can be treated symbolically: the lexicographic 
perturbation. 

We introduce the following notation that will be used for the reminder of the paper. 

Definition 5.1. The vector e := (e, e^, e^, . . . , e")^ is called the lexicographic perturbation vector 
and is a function of a positive real number e. 

We denote with S"^ the affine subspace S perturbed by e 

S' ■.= S + e = {y£W\y = Q6 + q + e,9eR'^} . (5.20) 

Theorem 5.2. Let M be a sufficient matrix and S be the affine subspace {Q6 + q\9 G M } for a 
given matrix Q and a vector q. There exists a 6 > such that S'^ := S + e lies in general position 
for each e G (0, 5). 
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Figure 2: An affine subspace 5 that is not in general position. (See Example 2) 

Remark 5.3. In the remainder of the paper we will use the standard expression "property A holds 
for all sufficiently small e > 0" rather than the more cumbersome "there exists 5 > such that 
property A holds for each e € (0, 6)" . Therefore the claim of Theorem 5.2 can be written as: S"" lies 
in general position for all sufficiently small e > 0. 

To prove Theorem 5.2, we need the following lemma. 

Lemma 5.4. Let Q G M"^'^, q G M" and S"" = {Q6 + q + e\ 6 ^ R"'}. For any complementary cone 
C, there exist finitely many e such that S^ intersects C but not int(C). 

Proof. Let -B be a complementary basis. We denote with hi the hyperplane hi := {y\(3iy = 0} for all 
i ^ B. Therefore C{B) n /ij is a facet of C{B) = {y £ W^\f3y > 0}. We will prove that for any subset 
J (^ B there are finitely many e such that the following condition holds 



^S%= C{B) nS' Q C{B) n (f^ hi) and S% 1 h for alH J 



(5.21) 



ieJ 



The statement of the lemma will then follow directly, since there are finitely many subsets of B. 

Let J be any nonempty subset of B. For any e for which (5.21) holds, J contains the indices of 
all inequalities of S"^ that are implicit equalities and it holds that 



Pj{QO + g + e) = for all ^ G TZ% 



(5.22) 



where n% = {e\Qe + q + e £ S%]. 

We can distinguish two cases. The case 1: there exists i G J such that (3iQ is a zero row vector. 
Since S^ is nonempty, /3j(g + e) = for the condition (5.22) to be valid. This non-trivial polynomial 
equation holds for at most n values of e. 

Now consider the case 2: l3jQ has no zero rows. We will prove that the matrix l3jQ does not 
have full row rank. Since 5^ is nonempty and has dimension d — rank(/3jQ) < d, the Chebyshev 
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center problem has an optimal value of zero, i.e. maxjt | fiiQO — t> —(3i{q + e) Vz} = 0. We consider 
its dual problem 

min {(3{q + e)Yy 

.X. -mfy = 

y > . 

From strong duality, there exists a non-zero optimal solution y* > such that {(iQ)'^y* = and 
{I3{q + e))-^y* = 0. We now claim that all indices of the strictly positive components of y* are 
contained in J. Assume that there is an index i ^ J with y* > 0. Then, for any 9 £ 7^^, 
f3i{Q0 + 9 + e) = 0, i.e. S'^ C hi, which contradicts the maximality condition in (5.21). Therefore, 
y} > and {(3Q) yj = 0, i.e. there exists a non-trivial combination of rows of PjQ. 

Let {vi, . . . , Vs} be a set of columns of Q such that /3jf i, . . ., f^jVs form a basis of the column 
space of (3jQ. Since PjQ does not have full row rank, s < | J|. If f3j{q + e),Pjvi, . . . ,(3jVs are 
linearly independent then for any G M'^ the equation (5.22) cannot hold. We claim these vectors 
are linearly dependent for finitely many e. 

First we consider the case s = | J| — 1. Then, these vectors are linearly dependent if and only if 

detipj{q + e),f3jvu...,i3jvs) = . (5.24) 

This condition is a polynomial equation in e and holds for finitely many e. Finally, if s < | J| — 1, 

we use the same argument by adding a proper number of vectors vi, . . . , ^ij|_i_s such that 

PjVi, . . . , (3jVs,vi, . . . , ^i j|_i_s are \J\ — 1 linearly independent vectors. D 

Proof of Theorem 5.2. We can assume without loss of generality that S does not lie in general 
position. For any complementary cone C exactly one of the following cases holds: 

1. S does not intersect C , 

2. 5 intersects the interior of C , 

3. S intersects the boundary of C and S n int(C) = , which can be differentiated into two 
subcases: 

(a) 3(5 > such that C n 5" = for aU e G (0, 5) 

(b) For all (5 > there exists an e G (0, 5) such that C n 5' / . 

For each of these cases we need to prove that there exists 5 > Q such that either S^ intersects int(C) 
for all e G (0,(5) or C n S"^ = for all e G (0,(5). Clearly, this condition holds for cases 1, 2 and 3a. 
Therefore, it suffices to prove that for case 3b there exists a (5 > such that S*^ intersects int(C) for 
alleG (0,(5). 

Let C be any complementary cone that satisfies condition 3b. From Lemma 5.4 and by assumption 
there exists a (5 > such that S^ intersects the interior of C and for any e G (0, 5) either int(C)n5'^ ^ 
or C n 5*^ = 0. More precisely, one can select any (5 > smaller than the smallest e > for which 5*^ 
intersects C but not int(C). Since S'' shifts continuously with e, there exists no e G (0,(5) such that 
C n 5^ = 0. D 

Example 3. Figure 5.1 shows two examples in which S does not lie in general position. In the first 
example (Figure 3(a)) the cones C({2,3}) = cone(/2,— Mi) and C({1,4}) = cone(Ii,— M2) contain 
adjacent critical domains, although they are not adjacent cones. In the second example (Figure 3(h)) 
two different critical domains coincide. In higher dimensions the critical domains can overlap in 
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(a) 




(b) 
Figure 3: Lexicographic perturbation of the affine subspace S 

several ways and therefore it is not evident how to choose an appropriate decomposition when this 
situation occurs. In both cases the affine subspace S can be artificially and symbolically shifted into 
general position through the use of lexicographic perturbation. 

5.2 Neighborhood computation in S"" 

Given a complementary basis B that is feasible in S^, the goal is to determine the adjacent 
critical domains to 5^. Since 5*^ lies in general position, it suffices to explore the adjacent bases of 
the basis B. Similarly to the non degenerate case, we first determine the facets of S^ at Line 5 of 
Algorithm 3 and then compute the adjacent critical domains that intersect with each facet. 

Let B he a complementary basis and consider S'^ = {y £ W^ \ I3y > 0, y = Q9 + q + e}, where 
P := A~Q . The hyperplane hi := {y \ /3iy = 0}, for some i £ B intersected with 5^ forms a facet of 
iS*^ if there exists a y* £ S'^ such that (3jy* > for all j G B\{i} and (3iy* = 0. Therefore hi fl 5^ is 
a facet of 3% if and only if 



r(e) 



max 

s.t. 



-(3jQ9 + t < (3jiq + e), yjGB\{i} 
-PiQe = I3i{q + e) 



(5.25) 



has a positive optimal value for all e > sufficiently small. 
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Algorithm 3 function neighbors^ (B): returns all bases whose critical domains are adjacent to 5^. 
Input: A complementary basis B, a sufficient matrix M and an affine subspace S^. 
Output: The set of complementary bases B, whose critical domains are adjacent to S^. 



1: 
2: 
3: 

4: 

5: 

6: 

7: 

8: 

9: 

10: 

11: 

12: 

13: 

14: 

15: 

16: 

17: 

18: 



B: 

(3: 
D 



A 



-1 

B 

A^ 



ixN 



■B ^-N G ^^ 

for each i ^ B do 

if isLexPosiii?;e(redundancy(i?,i)) is true then 
if D- > then 

add B' := B\{i} U {i} to B 
else 

for each j G B with D,^j < do 
B":=B\{i,j}U{i,j} 
if isLexPositive{ad^acency{B,B")) then 

add B" to B 
end if 
end for 
end if 
end if 
end for 
return B 



//where N = {!,..., 2n}\B 
//if Piy > is non-redundant 



/ /if dim{S% n S%„) = d-l 



This decision problem is no longer an LP, because the right hand side of the constraints depends 
on a polynomial in e and we want to know the behavior of t*{-) in the neighborhood of zero. In the 
next subsection we will propose an efficient method for determining if i*(-) is positive for sufficiently 
small e. 

If the hyperplane hi defines a facet of S'^, then we can distinguish the same three cases as 
discussed in Section 4.1. 



Diagonal pivot. If there is exactly one adjacent complementary basis B' , then S'^, is full-dimensional 
(by the general position of S"^) and is the unique adjacent critical domain to S"^ along S'^ D h^. 

Boundary of K{M). If there are no adjacent complementary cones to C{B) along hi (see Corol- 
lary 3.15), then there is no adjacent critical domain to 5^ along the facet S^ n hi. 

Exchange pivot. If there are adjacent bases B" with \B"f^B\ = 2 and dim{C{B)nC{B")) = d—1, 
then we must check for each such basis B" whether dim{S^„ n Sg) = d—1 (Line 11 of Algorithm 3). 
Assume B" := B\{i,j} U {i,j} where j € i?\{i}. As in the non-degenerate case, we formulate a 
decision problem similar to LP (4.19) to test the dimension of the intersection: 



t*(e) 



max 




t 






s.t. 


PkQe - 1 


> 


-Pk{q + e), 


yk G B\{i} 




(3^Qe 


= 


-m + e) 






(i'lQe - 1 


> 


-/?fe(9 + e), 


Vfc E B"\{j} 




P]Q0 


= 


-I3]{q + e), 





(5.26) 
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where /3 := A^ and /?" := A^„. The two critical domains are adjacent, i.e. dim(5|j n5'g„) = d- 
in S^ for all e > sufficiently small if and only if t*{e) > for all sufficiently small e > 0. 



5.2.1 Symbolic computation of the parametric Chebyshev center problem 

As seen in the previous subsection, the goal is to decide whether the optimal value t*(e) of (5.25) 
(and of (5.26)) is positive for all sufficiently small e > 0. We call this decision problem a parametric 
Chebyshev center problem. Here we introduce a method that can compute exactly the behavior of 
t*(-) for sufficiently small positive e; the proposed approach is summarized as Algorithm 4. The 
procedure is explained only for (5.25), since the same method can be easily applied for (5.26). The 
goal is to transform (5.25) into a multi-objective LP that can then be solved with any LP-solver. 
Recall the parametric LP (5.25): 

t*(e) := max t 

s.t. -/3jQe + t < Pj{q + e) , for all jeB\{i} (5.27) 

-(3iQ0 = A(g + e) . 

For a fixed value of e, this is a linear program and its dual is: 

t*{e)= min (/3(g + e))^y 

s.t. -iPQfy = 

Ej^^yj = 1 (5.28) 

Vj > for ah j G B\{i} 

Vi free . 

Let us denote its feasible region by F{i), that is, 

F{i) = {ye M"| - {(3 Qfy = 0, ^ y, = 1, yj > for ah j G B\{i}}. (5.29) 

In order to solve (5.28) symbolically we introduce the following standard notion. 

Definition 5.5 (Lexico-positive). A vector a G M* is lexico-positive (denoted by a ;^ 0), if a 7^ 
and the first non-zero component of a is strictly positive. Given two vectors x and y G M*, we write 
X y y if and only if x — y y 0. A matrix is called lexico-positive if all its rows are lexico-positive. 
If 5 = {s'lig/ is a set of vectors, then s^ is the lexico minimum of S if and only if s* >z s^ for each 

iel. 

The following theorem demonstrates that minimizing the polynomial cost function of (5.28) is 
equivalent to computing the lexicographic minimum of a vector. 

Theorem 5.6. If y is a feasible vector of the dual problem (5.28), then the two statements below 
are equivalent: 

1. 3(5 > such that {(3{q + e))'^y > for all e G (0, 5) , 

2. {l3[qI\fyyQ . 

Proof. The statement follows from the equality {[3{q + e)^ y = (1, e, e^, . . . , e")(/3[g /])^y, which 
holds for all y and for all e. For every polynomial p{x) = (1, j;,rE^, . . . ,x")(poiPi)P2> • • • ,Pn)'^ the 
following holds: there exists a, 6 > such that p{x) > for all x G (0, 6) if and only if the first 
non-zero coefficient of {po,Pi,P2, ■ ■ ■ ,Pn) is positive, i.e. {po,Pi,P2i ■ ■ ■ ,Pn) is lexico positive. D 
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We can now consider an equivalent problem that we call a lexicographic linear program (lexLP): 

redundancy(S,^):P*^= ^'^"^^ ^1^ I]V V^ ^ ^^^^ (5.3O) 

The cost of this optimization problem is vector valued and the operator lexmin means to compute 
the lexicographic minimum vector {f3[q, I])^y over all feasible decision variables y. We will denote 
this particular lexLP, which tests the redundancy of the i-th inequality in S'|j, as the function 
redundancy(i?, i) . 

Theorem 5.7. Ift*{-) is the optimal value of (5.27) as a function of e and T* is the optimal value 
of (5.30), then the following holds: 

36 >0 such that t*{e) > for all e € (0,5) ^^ T* y . (5.31) 

Proof. The statement is a direct consequence of Theorem 5.6. D 

Note that as is the case for linear programs, the restrictions and the objective function of (5.30) 
are linear, although the objective returns a vector instead of a scalar. We say that T* is the optimal 
value of (5.30). If the vector PiQ is non-zero then the feasibility region is bounded and therefore 
the optimal value is always attained if the problem is feasible. 

For our purposes, it is not necessary to compute the entire vector T*, but only a sufficient number 
of its elements in order to determine if it is lexico-positive or not. To this end, the goal is to find 
the first non-zero component of T* and therefore the lex min problem (5.30) can be treated as a 
multi-objective LP in the following way. First (say at step 0) we solve the LP: 

s.t. y e F{t), ^ ' 

where cq := /?(?. If Tq 7^ then we can conclude that the optimal value T* of the problem (5.30) is 
lexico-positive or lexico- negative from the sign of Tq. Otherwise, if To does equal zero, then we must 
consider the next objective function c\ := (i\ and minimise it while maintaining Tq = 0, and so on. 
If To = Ti = • • • = Tr-\ = 0, then at the step r we solve: 



min 


cjy 






s.t. 


y 


G 


m 




T 


= 


0, 



(5.33) 
A; = 0,...,r-1, 

where cq = P q and c^ = /3fc for A; = 1, . . . , n. If T^ 7^ is the first non-zero value of T* = (Tq, . . . , T„) 
then T* of (5.30) is lexico-positive if T^ > and lexico-negative otherwise. The resulting procedure 
is Algorithm 4, where the feasible region of the LP (5.33) is denoted by Fr. 
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Algorithm 4 Function isLexPositive{lexLP) 



Input: A lex linear program lexLP as redundancy(-, •) (5.30) or adjacency(-, •) (5.34). 
Output: Answer about lexico-positiveness of the optimal value T* of (5.30) or (5.34) respectively. 



1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 
9: 
10: 
11: 



Let Co, ci, . . . , c„ be the objective functions of lexLP and Fq be the feasibility region of lexLP 
for r = to n do 

t*. := minjc^y I U S Fr} //by solving the LP (5.33) 

if t; > then 

return T* is lexico-positive. 
else if t* < then 

return T* is lexico-negative. 
else 

F,+i := Frn{y\ cly = 0} 
end if 
end for 



Remark 5.8. Note that the lexLP adjacency(i?, i) always has a non-zero optimal value T* because 
zero is not feasible in (5.33) and the optimal solution y* of the last LP (5.33), with r = n, must 
be optimal also for all previous LPs. Since /? has full rank we must have that (3y* is non-zero and 
therefore there must be at least one component of T* that is non-zero. 

The parametric LP (5.26) that determines whether two critical domains are adjacent can also be 
solved using the same procedure. LP (5.26) can be rewritten as a lex min LP as follows: 

T* := lexmin {f3[q, I]f y + {(3"[q, I]fx 

s.t. _(/3Q)Ty_(^//g)T^ = 

'l2keB\{i} Vk + Z]fc6B"\{i} Xk = I 

yk > 0, keB\{i} 

Xk > 0, kGB"\{j} 
yi,Xj free. 

(5.34) 
We will call this lexLP adjacency(i3, B") because it tests the adjacency of Sb and 5^". Recall that 
one of the two equalities in (5.26) can be removed because one is redundant and therefore one of 
yi and Xj is also redundant. Hence, (5.34) has the same structure as (5.30) and can be solved as a 
multi-objective LP as explained above (see Algorithm 4). 



adjacency(i?,S'' 



5.3 Post-processing of the graph of critical domains Q"" 

The previous section introduced a computational method for computing the graph of critical do- 
mains G^ relative to the lex-perturbed space S^ for all e > sufficiently small. The goal in this section 
is to recover the graph G = (V, E) of critical domains relative to the original space 5 according to 
Definition 3.23. 

The following theorems will show that one can construct from Q'' the graph Q of critical domains 
relative to the unperturbed space S. 

Theorem 5.9. Let M he a sufficient matrix, S an affine subspace and consider the graph of critical 
domains G^ = {V", E^) relative to the lexicographically perturbed space S^ for sufficiently small e > 0. 
Let V be the set of all bases B in V^ with full- dimensional critical domains Sb- Then the following 
statements hold. 
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1. For each B G V" , the complementary cone C{B) intersects S and therefore Sb is nonempty. 

2. For any two distinct bases Bi and B2 in V , Sbi and Sb2 have disjoint relative interiors. 

3. The set of all full- dimensional critical domains Sb for B ^V covers Sf, i.e., 

\J SB = Sf = K{M) n s , 
Bev 

and forms a polyhedral decomposition of Sf. 

Proof. First, note that complementary cones are closed and so the first statement follows directly. 

To prove the second, note that for all e > sufficiently small (say e < (5) S^ lies in general 
position. From Proposition 3.25 the critical domains defined by S'' are disjoint in their interiors for 
any positive e < 6. The second statement follows from the fact that for any basis B, S^ D C{B) 
changes continuously in e, for e < 5. 

The third statement is proven in two steps. First we prove that the critical domains whose bases 
are in V^ define a covering of S" n K{M). Let q G K{M) n S, since e G K{M) and K{M) is a 
convex cone, there exists a basis B such that g + e G C{B) for all e > sufficiently small. Therefore 
q is in C{B) and hence in Sb because complementary cones are closed and thus Ufigy^ Sb = Sj = 
K[M) n S. Since critical domains are closed and Sj is a convex polyhedron, the full-dimensional 
critical domains define a covering oi Sj. D 

The above theorem demonstrates that the bases in V" have nonempty critical domains in the 
unperturbed space S. Moreover, there exists a subset V whose critical domains form a polyhedral 
decomposition of Sf according to Definition 3.23. The next theorem discusses how adjacency in Q'' 
relates to adjacency in Q. 

Theorem 5.10. Let M he a sufficient matrix and Bi,B2 be two complementary bases in V C V^, 
i.e. Sbi cmd Sb2 both have dimension d. If Sbi and Sb2 0,1"^ adjacent, then there exists a path 
B^, B^, . . . ,B^ in Q'^ = {V", E'^) from Bi to B2 with the following property: Sj^i intersects Sbj HSbj 
with dimension d — 1, i.e. dim(S'ji n Sbi H SB2) = d — 1 for all i = 1, . . . , r. 

Proof. If S^ and S"^ are adjacent, {Bi,B2) is clearly the desidered path. We assume they are 
not adjacent. We choose a ^ G relmt(SBi H SB2) which is not contained in any critical domain 
or in any face of dimension d — 2 or less, and let ^ G M"^ be the parameter with q = q + QO. 
We look now (for a moment) at the parameter space and at the critical regions. The hyperplane 
f := {9\ a 9 = b} contains the intersecion of the two critical regions, i.e. / ^ [T^Bi H T^Bj) ^^'^ 
consider its perpendicular (normal?) line 6{t) = 6 + ta. The image of 6{t) is q{t) = q + tQa in the 
original space S, respectively q'^{t) = q + e -\- tQa in the perturbed space S^. 

We know that for each e > sufficiently small every critical domain becomes either full- 
dimensional or empty. The full-dimensional ones vary continuously in function with e. Consider 
a segment [9*^(^1), 9*^(^2)] of the line {q'^it) \ t G M} such that it intersects either S'^ and S'^ for all 
e > sufficiently small. Because of the continuity no critical domain, which has dimension smaller 
than d — 1 in the original space S intersects this segment for all e > sufficiently small. Similarly, 
for any B & V^ no face of S'^ of dimension smaller than d — 1 intersects [(?^(ii), 9*^(^2)] for all e > 
sufficiently small. The desidered path is given by the critical domains which decompose the line 
segment between 5^^ and Sb2- D 

Note that the last condition in the above theorem, along with Proposition 3.21, implies that 
dim{Sj^i) = d — 1 for i = 2, ... ,r — 1 and therefore Theorems 5.9 and 5.10 imply the following 
corollary. 
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Corollary 5.11. Let M he a sufficient matrix, S an affine subspace and let Q^ = (V^,E^) be the 
graph of critical domains relative to 5"^. Then, the graph of critical domains Q = {V,E) relative to 

5 is related to Q'' as follows: 

1. y c y^ 

2. For every basis B ^V , Sb has dimension d. 

3. For each pair of bases Bi and B2, the critical domains Sbi o-nd 5^2 o-tc adjacent if and only 
if there exists a path {Bi, . . . ,Br) in Q^ with Bi = Bi, Br = B2 and dim{Sk) = d — \ for 
k = 2,...,r-l (or{Bi,B2) £ E' ). 

The above corollary provides a simple procedure for computing a critical region graph Q = {V, E) 
relative to the unperturbed affine set S from the perturbed one Q''. We begin from the perturbed 
critical region graph Q = {V,E) := G^ and remove each node B from Q that has a critical domain 
Sb which is not full-dimensional and add all new edges (i?i,i?2) to E satisfying the statement 3 of 
Corollary 5.11. From Theorem 5.9, the critical domains of the nodes of the resulting graph will form 
the desired polyhedral covering of the Sj. Theorem 5.10 states that the resulting graph contains 
edges for all adjacent bases, but may be overconnected since some of the critical domains Sb of 
removed bases B may have had a dimension less than d — 1. It remains, therefore, to test each 
edge in order to determine if the connected bases are in fact adjacent in the unperturbed space. As 
discussed previously, both operations for testing full-dimensionality and adjacency can be posed as 
linear programs. 

Remark 5.12. Note that much of the computation required to test for full-dimensionality of the 
critical domains for the unperturbed affine set has already been done while building the perturbed 
graph. Specifically, one can determine if a region is full-dimensional by examining the first component 
To of the optimizer of LP (5.30). 

6 Complexity of the algorithm 

In this section we will discuss the complexity of the proposed algorithm, which enumerates all full- 
dimensional critical domains relative to the lexicographically perturbed affine subspace 5*^. The 
well-known example by Murty (see [21] or see Chapter 6 in [22]), which was used to prove the non- 
polynomiality of the Lemke and the principal pivoting methods, can be easily seen to demonstrate 
that the number of critical domains of a pLCP with an affine subspace S of dimension 1 and P- 
matrix M is exponential in n. Since the complexity of the graph search (Algorithm 1) is a polynomial 
function of the number of critical domains, no algorithm for pLCP is polynomial in n. It is, however, 
possible to bound the number of operations required to explore the neighborhood of each critical 
domain, i.e. the complexity of Algorithms 2 and 3. Since each critical domain will be explored 
exactly once, we can say that the algorithm is output sensitive in that its complexity is a polynomial 
function of the number of full-dimensional critical domains and the size of input, provided that a 
polynomial-time algorithm for linear programming is used. 

We first consider the general position case and study the complexity of Algorithm 2. Assume that 
M is of order n, the affine subspace S is of dimension d and let B he a complementary basis with a 
nonempty critical domain Sb- The main computations of the function neighbors{B) (Algorithm 2) 
are: 

• Redundancy checking at Line 5 

(Solve LP (5.28) with n variables and d + 1 constraints.) 
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• Checking adjacency for the case of an exchange pivot at Line 11 
(Solve LP (4.19) with 2n variables and d+1 constraints.) 

We denote the time necessary to solve an LP in standard form by Tip[var, eq), where var denotes the 
number of (nonnegative) variables and eq is the number of equality constraints. The time necessary 
to explore a critical domain can then be bounded as follows. 

Theorem 6.1. Let M G M"^" be a sufficient matrix and assume that the affine subspace S lies 
in general position. For each complementary basis B with a nonempty critical domain Sb the time 
necessary to explore the neighborhood of Sb is bounded by: 

2 

n — Ti 
nTLp{n,d+l) + ^^—TLpi2n,d + l). (6.35) 

Proof. Redundancy checking requires the solution of LP (5.28) once for each of the n inequalities 
of Sb, which takes nTLp{n,d + 1) time. Adjacency checking by solving the LP (4.19) is necessary 
only in the case that the considered adjacent basis B" differs by two elements from the basis B and 

2 2 

since there are at most " "" such bases, the second term " "" Tip{2n, d+1) follows. D 

If S does not lie in general position, then the lexLPs (5.30) and (5.34) are solved instead of 
LPs (4.18) and (4.19). Each lexLP can be solved as a sequence of at most n + 1 LPs with the same 
variables and constraints (see Algorithm 4), which leads to the following complexity bound. 

Theorem 6.2. 7/ M € M"^" is a sufficient matrix, then for each complementary basis B with 
nonempty critical domain S'^ the time necessary to explore the neighborhood of S*^ can be bounded 
by 

(n^ + n) TLp{n, d + I) + !L^ TLp{2n, d+1). (6.36) 

The above theorems bound the complexity of "exploring" a basis of the output in Algorithm 1 
(Line 5). The condition at Line 6, which can be verified in time bounded by the logarithm of the size 
of the output, ensures that each output basis is explored exactly once. As a result, the complexity 
of the algorithm grows linearly with the size of the output and so is output sensitive. 

7 Example 

In this section we present a simple illustrative example that arises from control theory. Consider the 
following discrete time constrained linear time-invariant system: 

where x G M^ is the system state, x~^ is the successor state and u € M is the system input. A 
common method of control for this class of systems is Model Predictive Control, in which we solve 
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Figure 4: Polyhedral partition (left) and optimal input (right) for the control problem (7.37) 
at each point in time the following finite horizon optimal control problem: 



r{e) =min Y, iWQ^kWl + WRukWl) + WQ/^nI 

k=0 

subject to 



Xk+l 



I D^'^+Us'"'^ 



||2;fc||oo<5 , ||tifc-i||oo < 1 , V/c G {!,... ,iV} 



(7.37) 



where 9 is the current state of the system, the prediction horizon A^ is 5 and the weighting matrices 
Q, R and Qj are the identity. For high-speed systems, such as electric power converters, the goal is to 
solve the above quadratic program as rapidly as possible, in some cases at rates exceeding hundreds 
of kilohertz (e.g. [3]). By computing the optimizer offline as an explicit piecewise-affine function of 
the state 9, these speeds can be achieved [25, 14, 4]. The above parametric quadratic program is 
easily converted to a pLCP with a positive semi-definite matrix M [22], which was then solved using 
the proposed algorithm. The resulting polyhedral partition and mapping from the parameter 9 to 
the optimizer uo is shown in Figure 4. 

8 Conclusion 

In this paper an algorithm to enumerate all feasible bases of the parametric LCP defined by a 
sufficient matrix M and a lexicographically perturbed affine subspace S was proposed. It has been 
shown that the perturbed parametric LCP can be solved in a time linearly bounded by the size of the 
output and moreover, this output can be efficiently post-processed in order to generate a polyhedral 
decomposition for the unperturbed original affine subspace 5. 

One feature of the algorithm which is not ideal is the space requirement. Namely, the proposed 
algorithm must store all discovered feasible bases in the memory because it relies on the standard 
graph search technique. A great improvement can be made if we could apply the reverse search 
technique [1] which is essentially memory free. For this, it is necessary for the underlying graph 
to be oriented properly with exactly one sink. Somewhat similar to the present work, the paper 



24 



[9] proposed an algorithm to compute a polyhedral complex known as the Grobner fan which was 
shown to have such a "reverse search property." Finding such an orientation for the graph of critical 
domains is an excellent subject of the future research. 
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A Useful properties of matrix classes 

This section gives an overview of some matrix classes with important properties for linear comple- 
mentarity problems. The reader is referred to [6] for a thorough survey. 



P-Matrices 

Definition A.l. The matrix M € M"'^'^ is a P-matrix if and only if all principal minors of M are 
strictly positive. 

This class characterizes the matrices M for which the corresponding LCP always has a unique 
solution. 

Theorem A. 2. The following statements are equivalent: 

1. MgP, 

2. The LCP defined by the matrix M has a unique solution for all right hand side vectors q € M", 
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3. M does not reverse the sign of any nonzero vectors, i.e. 

[zi{Mz)i < for all i] ^ [z = 0] . 

Recall that M e M"^" is a positive definite matrix, denoted with PD if for all j; G M" it holds 
that X Mx > 0. It is then easy to see from the above theorem that positive definite matrices belong 
to the class P. 

Po-matrices 

Definition A. 3. The matrix M G M"^" is a Po-matrix if and only if all principal minors of M are 
non-negative. 

Analogously to the positive-definite case above, positive-semidefinite matrices (PSD) are clearly 
in Pq. The following theorem gives properties of PSD matrices relevant to the solution of pLCPs. 

Theorem A. 4. Let M G M"^" be a matrix and I be the identity matrix of the same order. The 
following statements are equivalent: 

1. Me Po, 

2. For each vector x ^ there exists an index k such that Zfc / and Zk{Mx)k > 0, 

3. (M + el) is aP matrix for all e> 0. 

Semimonotone matrices 

Definition A. 5. A matrix M G M"^" is called sem,im,onotone if the following holds: 

for all X > 0, rr / ^ [xk > and {Mx)k > for some k] . (1.38) 

The class of such matrices is denoted by £"0 and by Theorem A. 4, every Po-matrix is semimono- 
tone. 

Definition A. 6. Let M G M"^*^ be a semimonotone matrix. If for all index subsets a Q {1, . . . , n} 
with det{Maa) 7^ the principal pivot transform, of M with respect to a 



M' 



Ms^aM-^ Mr--i - Ms,^M-^M^s._ 



is semimonotone, then M is called fully semimonotone. The class of such matrices is denoted with 
Eq and M is said to be an E'q -matrix. 

Qo-Matrices 

Definition A. 7. An LCP defined by the matrix M and right hand side vector q is called weakly 
feasible if there exist positive vectors z and w such that w — Mz = q and feasible if z'w = also 
holds. The class of matrices M for which the LCP is feasible whenever it is weakly feasible, is 
denoted by Qq. 

Since P-matrices are feasible for each vector g, they are also Qo-matrices. 

Theorem A. 8. Let M G M"^" and L be the identity matrix of same order. The following statements 
are equivalent: 
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1. M^ Qo, 

2. The complementary range K{M) is convex, 

3. K{M) = cone ([I -M]) 

The implications of convexity of the complementary range are discussed in the next section. 

Sufficient Matrices 

Definition A. 9. A square matrix M is called column sufficient if it satisfies the implication: 

h(Mz)i < for alH] =^ [^^(Mz)^ = for alH] . (1.39) 

The matrix M is called row sufficient if its transpose is column sufficient. If M is both column and 
row sufficient, then it is said to be sufficient. 

Theorem A. 10. If M is row sufficient matrix, then 

1. M G Po, 

2. M e Qo- 

From the theorem above we have that every column sufficient matrix also belongs to Pq- Below we 
state a characterisation of column sufficient matrices, which has an important implication regarding 
the structure of the resulting complementary cones. 

Theorem A.ll. Given a matrix M G R"^", the following statements are equivalent: 

1. M is column sufficient, 

2. For each vector q € M" the following holds: if {z^,w^), {z'^,w'^) are two solutions of the LCP 
defined by the matrix M and the vector q, then {z^Y'w'^ = {z'^)'^w^ = 0. 
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